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ABSTRACT 

A 2D model representing the dynamical interaction of dust and gas in a planetary channel is explored. The two 
components are treated as interpenetrating fluids in which the gas is treated as a Boussinesq fluid while the dust is 
treated as pressureless. The only coupling between both fluid states is kinematic drag. The channel gas experiences a 
temperature gradient in the spanwise direction and it is adverse the constant force of gravity. The latter effects only 
the gas and not the dust component which is considered to free float in the fluid. The channel is also considered on 
an f-plane so that the background vorticity gradient can cause any emerging vortex structure to drift like a Rossby 
wave. A linear theory analysis is explored and a nonlinear amplitude theory is developed for disturbances of this 
arrangement. It is found that the presence of the dust can help generate and shape emerging convection patterns 
and dynamics in the gas so long as the state of the gas exceeds a suitably defined Rayleigh number appropriate for 
describing drag effects. In the linear stage the dust particles collect quickly onto sites in the gas where the vorticity is 
minimal, i.e. where the disturbance vorticity is anticylonic which is consistent with previous studies. The nonlinear 
theory shows that, in turn, the local enhancement of dust concentration in the gas effects the vigor of the emerging 
convective roll by modifying the effective local Rayleigh number of the fluid. It is also found that without the f-plane 
approximation built into the model the dynamics there is an algebraic runaway caused by unrestrained growth in the 
dust concentration. The background vorticity gradient forces the convective roll to drift like a Rossby wave and this 
causes the dust concentration enhancements to not runaway. 



1 INTRODUCTION AND SUMMARY OF RESULTS 

There has been growing interest in recent years in the possibility that persis- 
tent vortex structures in protoplanetary discs may be the sites where plan- 
etesimals are formed (Barge & Sommeria, 1995, Tanga et al, 1996, Bracco 
et al, 1999, Barranco & Marcus, 2000, Barranco & Marcus, 2006). These 
investigators have put forth the scenario in which a disc, laden with dust, 
supports some type of long-lasting vortex (driven by some unspecified exci- 
tation mechanism, except for Klahr & Bodenheimer, 2003 and Barranco & 
Marcus 2006) that manages to attract the dust through the combined action 
of gas-particle drag and Coriolis forcing (for example, Barranco & Mar- 
cus, 2000). Separately, it is also well-known that strongly shearing flows 
(Keplerian discs qualify under this classification) in rotating frames pref- 
erentially support the persistence of anti-cyclonic structures while almost 
entirely wiping out cyclonic vortices (Bracco, et al, 1999). In light of 
this, one of the interesting results of the aforementioned dust-vortex inter- 
action studies is that anti-cyclonic vortices also happen to be the sites onto 
which disc dust collects most rapidly (most notably, Bracco et at, 1999). 

In the investigations mentioned, the dynamics of these dust-disc sce- 
narios are treated as one-way: the dust passively responds to the gas flow 
without any back-reaction of the dust upon the gas. Given that the current 
hypothesis of protoplanetary discs is that their dust contributes no more 
than a few percent of the total mass density of the disc (Hayashi, 1983), it 
is reasonable to treat the dust as a collection of Lagrangian tracers which 
responds to the vortex induced gas flow and having no dynamical influence 
on the gas itself. 

However, it is interesting and instructive to turn the question around 
and examine what would result if the dust does dynamically effect the gas. 
One could envision a situation in which either the gas in the disc has been 



removed sufficiently or one is investigating the dynamics taking place near 
the disc midplane, where planetesimals are likely to be strongly concen- 
trated. 

As a thought experiment, suppose the gas in a model shearing envi- 
ronment develops, through some type of dynamical instability, a series of 
long-lived vortex structures. In addition to obviously effecting the dust tra- 
jectory in the way it is usually envisioned, it seems reasonable to suppose 
that the presence of dust would dynamically influence, because of their mu- 
tual dynamic coupling (e.g. kinematic drag), the manner in which emerging 
gas structures evolve. This all would be plausible under conditions in which 
the gas and dust have equal dynamical influence upon each other. Youdin 
and Goodman (2005) have conducted a similar study considering the linear 
instability emerging from interpenetratingstreams of gas and dust whose 
only interaction is through kinematic dragQ 

Presented here is a simple and idealized model for the way this type of 
interaction develops in a model confined environment. The physical model 
employed here is a flow in a rotating channel. The material is treated as two 
interpenetrating fluids of uniform density, one representing a pressureless 
dust "fluid" and the other an incompressible Boussinesq gas. The gas ther- 
modynamics is governed by thermal conduction and we suppose that there 
is a constant gradient of the gas temperature from one wall to the other (in 
the spanwise or "wall-normal" direction). As part of the idealization that 
goes into the model considered here, the gas component alone is subject to 
a constant force in the wall-normal direction and that the Boussinesq buoy- 
ancy effects are associated with this force in the usual way. We suppose, 
further, that the gas and dust are both subject to an external force which 
gives rise to a linear Couette shear of both fluids in the direction parallel to 
the channel walls. The coupling between the gas and dust is through a sim- 
ple kinematic (Darcy) drag prescription like, for example, in Tanga et al. 



1 A coherent vortex structure is said to be anti-cyclonic if the signature of 
the background vortex state in which it is embedded is opposite the sign of 
the vortex structure itself. 



2 There have been a number of studies of this sort in which the interpene- 
trating streams are coupled to each other through gravity. For a more recent 
study see Bertin & Cava (2006) and references therein. 
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(1995). There is no viscosity. The ingredients of this hypothetical scenario 
already predisposes the gas component susceptible to buoyant instability ala 
Rayleigh Benard. Finally, because the dust fluid is treated as being pressure- 
less, no a priori restrictions upon the compressibility of the dust is made. 

In Section 2 we show the development of Rayleigh-Darcy type of con- 
vection. In the formulation of the problem, the dust density decouples from 
the linear analysis of the developing convection roll. A simple relationship 
is demonstrated that further shows that a steady vortex roll at the corre- 
sponding critical wavenumber causes the dust density to grow algebraically 
when the roll is anti-cyclonic and, conversely, the density is depleted where 
the roll is cyclonic. 

In Section 3 a nonlinear asymptotic analysis is developed in the limit 
of large aspect ratio and under conditions of fixed thermal flux on the chan- 
nel walls and where the background vorticity is modeled as an f-plane with 
a weak gradient. In addition to this, the imposition of a weak external force 
will promote a steady flow exhibiting a weak amount of shear (here it will 
be Couette). The asymptotic analysis proceeds under the assumption that 
both (i) that the thermal time is much greater than the dust stopping time 
and, (ii) that the time scale derived from the geometric mean of the ther- 
mal and rotation times is also much greater than the stopping time of the 
dust. The two conditions translates to a situation in which the dust stop- 
ping time is actually much longer than the local rotation time. This scaling 
regime also implies that the dust velocity behaves as if it were irrotational 
at leading order. 

The model shows algebraic instability (in the dust component) unless 
some amount of background vorticity gradient is present in the flow (as pro- 
vided by the weak f-plane prescription). This gradient stabilizes the growth 
by promoting a Rossby wave drift of the the convective roll (vortex) pat- 
tern. Because the Rossby wave drift speed of the vortex pattern is different 
than the background velocity field the effect here is for the fluid pattern to 
drift past the dust component. As given regions of the dust are exposed to 
overall reductions of the total fluid vorticity (as measured with respect to 
the background vorticity arising from the channel's rotation) the dust be- 
gins to collect as expected from the behavior predicted by the linear theory. 
However because the pattern is not stationary the dust cannot locally collect 
ad infinitum because the low vorticity part of the vortex pattern will pass on 
through a given local dust region. In turn, the dust region will be exposed 
to patches of fluid passing by with enhanced vorticity and, consequently, 
will experience a reduction in the local dust concentration. It is found, how- 
ever, that secondary growing oscillations also appear in this formulation and 
these are wiped out when a certain amount of diffusion in the dust concen- 
tration is added to the model. It is also found that according to the model, 
places where the dust concentration is enhanced (depleted) corresponds to 
zones where the convective roll is slightly weakened (strengthened). We 
show that this can be explained by an effective modification of the local 
Rayleigh number due to the nonlinear rearrangement of the dust concentra- 
tion. 



2 TWO DIMENSIONAL FLOW EQUATIONS 

Consider a two fluid description of a dust laden simple Boussinesq fluid in 
a rotating channel. For the sake of generality we suppose that only the gas 
phase has certain properties which cause it to be dynamically influenced by 
an external constant force of magnitude g and with direction normal to the 
channel walls (viz. in y direction). In this toy model we additionally posit 
there to be a second force represented by the term which also acts in the 
y direction and its dependence is on y only: in other words, = (y). 
Unlike the other force g, the force affects both gas and dust. The fluid is 
dynamically coupled to the dust because of simple kinematic drag measured 
by the difference in velocities between the gas and dust phase. Finally, we 
assume that the thermodynamic transfer properties of the gas are modeled 
by simple thermal conduction. Thus, the equations of motion for the gas 
phase are summarized to be, 
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In which u is the gas velocity and u d is the dust velocity, P is the gas 
pressure, The dust density is p d while the gas density is simply p. The 
superscript (x) and (y) appearing in front of velocity expressions are re- 
spectively their x. and y components. The kinematic drag between phases 
is mitigated by the "stopping-time" t s and we assume it to be constant - 
see similar treatments by Cuzzi (1993) and Tanga et al, (1995). If the chan- 
nel rotates with rate Q in direction z then the Coriolis term / is given by 
2Q. Because the gas phase is assumed to be Boussinesq its thermodynam- 
ics is modeled by the usual expression of a fluid undergoing fluctuations 
under constant pressure conditions where C p is the specific heat at constant 
pressure. 

Since this is a two-dimensional model, it will be more convenient for 
us to consider the gas phase momentum equations in terms of a single equa- 
tion for the gas vorticity. Since the gas phase is divergence free except when 
coupled to g, we may consider the gas velocity as resulting from a stream- 
function tp such that 
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It follows that by defining the vorticity as the curl of the fluid velocity, 
u> = V X u, then we have 
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thereby eliminating direct reference to the gas phase's divergence free na- 
ture and reducing the number of equations from three to two. The Jacobian 
appearing in {6) is defined as 
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Unlike its gas counterpart, the dust phase is assumed to be compressible, 
and this is why we have retained all spatial dependences of p d in (6). Fur- 
thermore, the dust phase vorticity, u) d , is defined from the dust phase veloc- 
ity, u d , in the usual way, viz., 



du 



(y) 



du 



(x) 



^d 



dx 



dy 



(7) 



The dust phase is modeled as a pressureless, compressible fluid under 
the influence of particle drag, the force and the Coriolis effect, 
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Since the underlying mass density states are presumed to be constant in 
both fluids, it will prove beneficial to write departures of the dust density in 
terms of a fractional density parameter, 8, defined by, 

Pd = p d {l + S). (10) 



2.1 Steady State 

A thermal steady state is when the gas conducts in the y direction. This 
means that there is linear temperature gradient in y given by, f = T* — 
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fixy, in which T* is the background temperature scale and where the con- 
stant /3y represents, in a sense, the steady background heat flux of the en- 
vironment. Denoting u and u d as the steady velocity profiles, a solution is 
u = u d , along with, 
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And as previously mentioned, the densities of both fluid states, p and p d , is 
assumed to be constant. We consider steady conditions to have no x direc- 
tion dependences. Thus, the temperature steady state condition is 

,d 2 f 
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(12) 



which yields a steady flux state of the gas with a temperature profile f = 
T„ — /3y, where /3 is a constant related to the background flux state of the 
gas and where T» is the background temperature scale. 



2.2 Linear Theory 

For the sake transparency, we consider linear disturbances of the steady 
arrangement prescribed above when the force Tj, is zero. This means that 
perturbations ensue under static conditions: u = u d = 0. Temperature 
fluctuations are represented by 



T = T, - f3y + 9, 



(13) 



where 8 represents the temperature disturbance away from the static state. 
Except for buoyancy effects, the gas density is constant. Therefore the 
Bousinessq approximation says that, 



(14) 



where here the coefficient of thermal expansion at constant pressure is rep- 
resented by a. A prime on any density quantity is meant to designate its 
deviation from the steady uniform state, which will be designated with an 
overbar over the quantity. 

We suppose that all dynamical disturbances occur on a length scale 
of d. Consequently spatial scales are nondimensionalized on d and, further- 
more, they are written as 
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t = Cppd 2 1 K is a thermal time scale and, as such, we scale all temporal 
derivatives by it. Temperature disturbances are scaled by fid and is written 
in nondimensional terms as 9 = @(3d. Velocities scale with d/r which sets 
the scaling for both the stream function and vorticity (for both fluid phases), 
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in which tfi and w are the nondimensional stream function and vorticity, 
respectively. In this formalism the nondimensional linearized perturbation 
equations for the gas are, 
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Whilst for the dust phase, 



d_ 

at 



Pd t 

1 Quid 
Pr at 
i ac 
Pr at 



o. 

-(uj d - u) - fC, 
-C + fuj d . 



(15) 
(16) 
(17) 

(18) 
(19) 
(20) 



The variable C, defined by 
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represents the dilatation of the dust phase. The term u d represents the dust 
velocity scaled to d/r. The Rayleigh number, R, is defined to be 
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while an effective Prandtl number Pr is identified with 
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Finally, the Coriolis parameter is scaled to the stopping time: / = ft s . 

We note that in the dust momentum equation, out of which both H9\ 
and i20\ are formed, all appearances of the dust density and its fluctuation, 
p' d /p d , explicitly drop out, In practice this means that (TF) explicitly decou- 
ples in a normal mode analysis of the perturbation equation set )15M20t . 
This observation has an interesting consequence for the nonlinear theory in 
the upcoming section. 

The limited coupled set of equations I15M17I and U9M201 admit 
steady (i.e. d/dt = 0) solutions. The particle phase variables are simply 
related to the gas phase variables by, 
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Note that {25) says that the dust dilatation is related to the gas vorticity, 
which will have consequence shortly. The gas phase equations become, 
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along with )17t unaltered as it appears. Combining U7t with 124127) yields 
a single equation for ifi, 
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(28) is supplemented with boundary conditions. As for the thermal condi- 
tions: we explore two different types of requirements, namely that the (i) the 
temperatures are fixed at the boundaries or (ii) the thermal fluxes are fixed 
at the boundaries. The former condition is the more familiar of the two in 
which, in practice, it requires that the temperature fluctuations are zero at 
the boundaries. The second condition is, effectively, thermal insulation and 
it, in practice, amounts to setting 80 /dY = at the boundaries. 

Assuming, in general for all quantities, periodic solutions in the hori- 
zontal £ direction with corresponding horizontal wavenumber q, e.g. 

ifi = "ip(Y) expiq^ + ex., 

we find that for fixed-temperature boundary conditions the lowest normal 
mode solution in Y for ■tp is 



ip ~ exp iqi; cos ttY, 



(29) 



where the system eigenvalue, which is an effective Rayleigh number defined 

by, 
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Figure 1. Critical Rayleigh number as a function of scaled horizontal wave- 
length q. (a) Conducting boundaries and (b)insulating (fixed-flux) bound- 
aries. 



There is a minimum value of R e ff at a critical wavenumber q c required 
for such marginal solutions to exist. We denote the minimum value of the 
effective Rayleigh number as R%* and we find that for fixed-temperature 
boundary conditions q c = tt and, R% = 4tt 2 «i 40 - which is the 
number expected for Rayleigh-Darcy convection (Horton & Rogers, 1945). 

The solution for perfectly insulating boundaries resists the straightfor- 
ward form characterizing the fixed temperature result. It proves to be easier 
to seek numerical solutions here instead. Figure[T]b) displays the Rayleigh 
number as a function of horizontal wavenumber for these marginal condi- 



tions. It turns out that K; 



12 and it occurs at q = 0. The result that 



the critical Rayleigh number occurs at infinite horizontal scales, though cu- 
rious, is a well known result in standard Rayleigh-Benard problems and 
its implications have been studied by many investigators (e.g. Chapman & 
Proctor, 1980). 

Irrespective of the thermal boundary condition assumed the marginal 
roll state can cause the dust density to grow algebraically. To see this con- 
sider \2Q\ with the dust dilatation term re-expressed in terms of the rela- 
tionship |25j. It is found that the time rate of change of the dust density is 
expressed by 



dt \ p d 
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(32) 



In other words, the dust density changes linearly in time (since there is 
no time dependence in w) and, additionally, this change is fastest at places 
where the vortex amplitude is extremal. Given the relationship in i32\ it is 
apparent that dust density grows fastest near anti-cyclonic vortex perturba- 
tions : this is clear since the RHS of the expression is positive only if the 
product — fu) is positive and this can only happen if the signs of / and u) 
are opposite. 

To be more accurate in this description of growing dust concentration 
it should be stated that if the sign of the growing vortex perturbation is 
negative it only means that this part of the fluid state starts appearing to have 
lower vorticity than the original global state. The result of this linear theory, 
taken in this light, is not too surprising because earlier work on the behavior 
of Lagrangian tracer particles in fully developed isotropic turbulence flow 
(Squires and Eaton, 1991) show that particles concentrate in those parts of 
the flow which are regions of lower vorticity and high strain rate. 



3 AN ASYMPTOTIC NONLINEAR MODEL 

In order to understand how to handle the curious situation involving the 
algebraic growth of the particle density perturbation even though the gas 
phase is neither growing nor decaying we asymptotically analyze the gov- 
erning equations under certain extreme conditions. In particular, we will 
explore the nonlinear development of disturbances for fixed thermal flux 
boundary conditions. Furthermore we will assume that the Prandtl number 
is much larger than the scaled Coriolis parameter, / which will also be much 
greater than one, viz, , 






(d) Maximum Absolute Amplitude, Time Series 





Figure 2. Representative profiles for Model A: /3 = 61 = and = 5.30 
and <fi = 10. Plotted are the lowest order temperature amplitude F, the 
lowest order vorticity amplitude F x , followed by the dust concentration, 
A. Panel (a), at time s ~ 0.21. Panel (b), time s ~ 1.25. Panel (c), blow 
up time s ~ 1.87. The scaled dust concentration always seems to trace 
the gas vorticity F x even throughout the blow up time. Panel (d) is the 
time series for the absolute values of vorticity scale and dust concentration, 
respectively. 



This extreme ordering says physically that the dust stopping time, t s , is 
much longer than the local rotation time ~ Q ~ 1 . This scaling however also 
says that the dust stopping time is much less than the geometric mean of the 
thermal and rotation times. In other words, 



n 



and 



(it s > 1. 



For the sake of this asymptotic calculation we take Pr — » 00. 

Furthermore we will suppose that the horizontal lengths scales are 
much larger than the scales in direction of gravity and we will introduce the 
small scaling parameter e to represent this extreme state. Corresponding to 
this situation we introduce the scaled horizontal variable X = e£, where X 
is an order one quantity. Additionally we introduce an order 1 scaled time 
variable T such that T = t 4 t. Thus, we replace all derivative operators 
with 
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From now on we will appeal to the shorthand notation for derivative op- 
erations for compactness of notation and clarity of derivation. Correspond- 
ingly, we find that a distinguished limit exists when the stream function is 
O (e) the temperature fluctuation. Specifically, we say that 

■0 = e*, lo = efl, 

where ^ and Q are now order 1 quantities. 

We introduce the external force, T^, into the analysis in which its 
nondimensionalization is expressed through the relationship, 



where jF fe is nondimensional. We assume that the steady external force has 
the form 



F k = -b<j>Y, 



(34) 



Pr > / > 1. 



where b and <h are order 1 constants. 
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In order to achieve proper asymptotic balance we exploit the freedom 
we have in choosing the "largeness" of the Coriolis parameter by scaling it 
(ad hoc) as 
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where cj> is, as before, order 1. The bizarre choice behind this asymptotic 
ordering is motivated by introducing the dust disturbance term S into the 
subsequent analysis at the right stage. The j3Y term, which tacitly repre- 
sents an assumption of a weak f-plane effect (see Pedlosky), is introduced 
a priori in order to ensure stability of the nonlinear model. Though it is an 
artificial device, it has the effect of causing a Rossby wave type of drift of a 
vortex pattern. 

With respect to jilt , the resulting steady x-direction flow is weakly 
Couette. Denoting the steady flow by U and U d then we have 



U = 



U d = ebY + O (e 7 ) 



(36) 



From here on out, we consider disturbances of the dust and gas velocities 
on top and over this base state. 

With these scalings in place we find that the dust quantities require the 
following scalings in order for there to be nontrivial balance: the fractional 
dust density scales as e 2 and from here on out it is written as 
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We find that the dust dilatation is first nontrivial at order e e and it is written 

as, 



C = e 6 C 6 + O (e S ) 



(38) 



whereas the dust vorticity scales even smaller such that it is nontrivial at 
order e 11 ! In other words, 



Ud - 



(39) 



Because of this last observation we will, for all intents and purposes, treat 
the dust velocity as derived from a potential which similarly scales on order 
e 6 . In other words, u d = e 6 Vr. As such, each component of the dust 
velocity scales as 



e 6 v 6 
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Together with the scaling relationship (38} and general definition 12 1 1 it is 
clear that at these initial orders the dust velocity satisfies a Poisson type of 
equation, 

?2l 
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With these scalings and assumptions in place we consider the temporal 
development of disturbances when the Rayleigh number deviates from its 
critical value by an order t 2 amount written as, 



R^R + t 2 R 2 . 



(42) 



For this analysis we will consider the development in terms of fixed-flux 
boundary conditions which means that the critical Rayleigh number about 
which we will be expanding around will be 12. However, the following 
calculation self-consistently selects the proper value of i?,o- The governing 
equations (in the Pr — > oo limit) with these scalings inserted now reads 
for the gas phase, 
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e i d T e + e 2 J(e,^) + e 2 b 1 Yd x e = e 2 dx* + e 2 d 2 < + Ui$) 

Q = (D 2 + e 2 a|)*, (45) 

and for the particle phase they are, 

= eQ - e<t>C 6 + O (e 7 ) (46) 
= ~e 6 C 6 + e 6 cf>n d + O (e 12 ) (47) 
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The new Jacobian symbol, J is defined on the variables X and Y . In the 
next section we set out to solve, order by order, the set of equations 1431481 
with fixed flux thermal boundary conditions. 



3.1 Expansions 



The set of equations 143148) are solved as a perturbation series expansion in 
powers of e 2 . Thus, the following expansions are presumed for this purpose, 



e = eo + e 2 e 2 + € 4 e 4 + - 
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<44t to lowest order becomes simply, 
D 2 o = 0. 
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The solution of this equation subject to the condition that the flux remain 
fixed on the boundary is, 



e = f(x,T), 



(54) 



meaning to say that the lowest order temperature perturbation is indepen- 
dent of the vertical coordinate. The rest of the expansion procedure is stan- 
dard and we relegate its full exposition to Appendix lAl What results are 
two coupled evolution equations for the temperature disturbance / and the 
lowest order particle density 5 2 in IA231 and )A32t . These are reproduced 
below, 
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The positive constants re, fi, B, S are defined in Appendix [A] r 2 is an ef- 
fective departure from the critical Rayleigh number i?,o = 12 given by 
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The parameter r 2 can be interpreted as representing an effective departure 
from criticality which includes the fact that a constant (linear) shear is a 
stabilizing influence against the onset of rolls in the linear regime. 



3.2 Canonization and Identification of Effects 

It is more convenient to define a number of rescalings of (55} and (56} in 
order to transparently display the structure of the two coupled equations. 
One may define new space and time coordinates, £ and s along with new 
amplitude scalings F and A by, 



X^y^X 

1/2 



f(X, r) (") F (X, s), 6 2 {X, T) -» A( X , •). 



(57) 



Consequently, {55} and {56} rewritten in terms of the definitions (57} is 



7> 3 ) 

x>x 



(ra - A)F, 



-6i(F 2 ) x -/3F xxx , 
-fax,, 



(58) 
(59) 



where 



<PRo 



The model set {58} and (59} are now in a more transparent form and allows 
us to discuss the effects they contain. The role of the first two terms on the 
RHS of (58} are well known (see Chapman & Proctor, 1980) respectively, 
(i) to damp the temperature profile due to horizontal thermal dissipation 
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(a) Time Series: no dissipation (b) Time Series: with dissipation 




5 10 15 20 5 10 15 20 

Time: s Time: s 



(c)s=0.84 (d) s =1.722 (e) s =20.958 




-2 2 -2 2 -2 2 



space space space 

Figure 3. Representative profiles for Model B: b± = and = 5.30, 
<t> = 10, = —10. (a), maximum amplitude of the dust concentration 
and the scaled vorticity when there is dissipation (B = 0). The amplitudes 
oscillate with growing size on a long-time scale, (b), maximum amplitude 
of the dust concentration and the scaled vorticity with dissipation present, 
(B = 1). The oscillations are quelled while the overall average maximum 
amplitude of the structures are preserved. Panels (c-e) representing profiles 
of F, F x , A at various stages in the ensuing evolution. It should be observed 
that the dust concentration, A, eventually locks onto the temperature profile, 
F, at late times. 



and (ii) to control amplitude due to nonlinear advection of temperature. The 
r'2 expression in the third term on the RHS of 458t constitutes the usual 
source of amplitude growth because of positive departures from the critical 
Rayleigh number. The A term expresses the fact that as the particle concen- 
tration grows the local Rayleigh number of the fluid goes down somewhat. 
This is best seen by observing the definition of the Rayleigh number (22): 
if the dust density goes up then R goes down and vice versa. The other 
terms on the RHS of 458t represent, respectively, (iii) the nonlinear twist- 
ing of the temperature profile due to shear and (iv) an effective Rossby 
wave type of drift of an emerging profile. As discussed in the previous sec- 
tion, {59j represents the local enhancement of dust concentration due to the 
emerging convection roll's vorticity, which is, to lowest order, proportional 
to F x as demonstrated in AppendixlAl 

Note that the nonlinear term {F x ) x associated with the shear is 
known to promote subcritical transitions in Rayleigh-Benard convection 
(Gertsberg & Sivashinsky, 1981) and will do so here in this model as well 
but we will not explore this feature in this study. This term appears geneti- 
cally so long as there is some sort of symmetry-breaking effect in the system 
(e.g. see Depassier & Spiegel, 1982). 



3.3 Selected Solutions 

The reduced asymptotic equations 158159) are evolved according to the nu- 
merical spectral scheme outlined in Appendix |B1 and is similar to the tac- 
tics used in Umurhan & Regev (2004). The robustness of the scheme was 
checked against the analytic solutions of a similar evolution equation de- 
rived by Chapman and Proctor (1980). The tests demonstrated agreement 
between the numerical scheme and the analytical results to one part in 
10 -7 . The results reported here were generally computed on 256 Fourier 
modes and when convergence was doubted, the results were checked against 
simulations based on 512 and 1024 modes. Because there are fourth order 
derivatives in the linear terms of 1581 we note that, (a) no artificial viscos- 
ity was needed and, (b) small time steps on the order of dt = 10 — 4 were 
typically taken. All simulations started out with random white noise for F 




Space 



Figure 4. Position of the peaks of the vorticity profile for Model B. The 
maximum peak is labelled with a circle and all secondary peaks are denoted 
with x's. After the transient readjustment phase is over (near s ~ 2) the 
system goes from multiple peaks to only one. After this time, the path of 
the peak is traced out with a dashed line and it is evident that it propagates 
from left to right on this spacetime diagram. 



Temperature, s-18.9 




Figure 5. Representative contours for Model B when the growth has settled 
onto the steady state. For these contours e = 1. (a) Temperature fluctuation 
contours, ©. (b) Stream function contours, 

whose total amplitude did not exceed 10~ 2 . The dust concentration always 
started out with zero initial amplitude. 

We consider here four selected models to discuss; 

• Model A. No background vorticity gradient and no shear, 62 = /3 = 0, 
while r'2 = 5.3 and the dust coupling is set <f> = 10, 

• Model B. The background vorticity gradient is negative, /3 = —10, 
there is no shear 62 = 0, along with T2 = 5.3 and <f> = 10 and, 

• Model C Some amount of shear with a prograde sense, bi = —10.0 
along with the same values of T2,4> and f3 of Model B. 

It should be noted that in all models the background vorticity is getting 
smaller with increasing Y. Additionally, the results reported here required 
there to be some dissipation added to the equation describing the dust con- 
centration evolution or otherwise the solutions blow up. This is discussed 
in more detail for Model B. 



3.3.1 Model A 

This and other models like it, where /3 = 0, are unstable. The amplitude 
and dust concentration steadily runaway once the dust concentration locks 
into the standing vorticity profile that emerges. To see this more readily re- 
fer to Figure|2] At early times in the simulation the local dust concentration 
grows fastest wherever G = F x , which is proportional to the leading or- 
der vorticity, is greatest. This behavior is in line with the conclusions of 
the linear theory section. For T2 = 5.3 and <f> = , the fastest growing 
wavenumber is n = 2 and this readily shown in the linear theory analysis 
of (58). Despite the linear theory expectations however, the final roll profile 
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will go from being double humped to being single humped as was shown by 
Chapman and Proctor (1980). This is because the nonlinear stability anal- 
ysis shows that the system prefers to eventually settle onto a profile with a 
single roll. 

With (f> non-zero, the same sort of late time behavior is observed. 
The dust concentration also qualitatively resembles the vorticity distribu- 
tion and it, similarly to the linear theory, appeal's to grow/deplete with- 
out bound. Those places where the vorticity amplitude is greatest are those 
places where the dust concentration is greatly depleted. Consequently, those 
locations in the gas correspond to greatly enhanced values of the local 
Rayleigh number thereby causing the roll profile there to continue grow- 
ing, in what appears to be, without bound. Naturally, the validity of this 
predicted blowup should be questioned since unbounded growth means the 
asymptotic solution has ceased being valid. 



3.3.2 Model B: Vorticity Gradient with No Shear 

Introduction of stabilizes the runaway growth seen in Model A. The am- 
plitude of the emerging convection roll settles onto the vicinity of a given 
value though on oscillation sets in. On a very long time scale this oscil- 
lation, which shows up in both the roll amplitude and dust concentration 
amplitude, exhibits a growth in amplitude and it eventually blows up as is 
evidenced in Figure[3^. This happens for all models in which <f> is not zero 
and the secular growth time scale is shorter as <j> becomes larger. In order 
to stabilize the long time behavior a dissipation term is introduced into the 
dust evolution equation so that, and from here on out, instead of i59\ the 
following equation is evolved, 

A s = -4>F X +t3A xx , (60) 

B is the scale of the dissipation and in the simulations it is taken to be 
unity. With the simulation run with £5 = 1, the resulting maximum absolute 
amplitude is shown in Figure [3p. The oscillations observed for B = are 
removed and the overall average amplitude of the underlying structures are 
preserved. 

The time evolution of these simulations begin with an initial growth 
phase. For T2 = 5.3, the linear theory predicts that the fastest growing 
Fourier mode is k = 2 and, as such, in the early stages the temperature 
profile is double humped as in Figure [3p. As Chapman and Proctor (1980) 
demonstrated, this type of profile is nonlinearly unstable and eventually the 
system undergoes a phase of transient readjustment (around s = 1.5) in 
which the temperature and vorticity go from being a two-humped profile 
down to a one-humped profile. Figure [5J; shows the state of this profile 
readjustment near the time s ~ 1.7. 

Additionally, Figure [4] demonstrates how the profile drifts with in- 
creasing x as a function of time. In this steady drifting state, the profiles are 
single humped and, as in Figure |5J, after the post readjustment phase the 
dust concentration profile locks onto the temperture profile tightly. The rea- 
son that the dust concentration does not grow without bound when the roll 
pattern drifts is simple. In Model A, because = 0, the temperature and 
vorticity profile that develops does so in place. Because the dust concentra- 
tion responds to local variations of the vorticity, places in the growing roll 
profile where the vorticity is most reduced are also where the rate of dust 
collection is greatest. Without a suitable saturation mechanism the greatly 
enhanced dust concentration begins to modify the local Rayleigh number of 
the fluid. This, in turn, causes the roll amplitude to grow with ferocity. With 
a larger amplitude in the vorticity will come an ever increased concentra- 
tion of dust. In this way the cycle feeds on itself and it runs away. On the 
other hand, when the roll profile drifts, local enhancements of dust do not 
grow without bound because the low vorticity patch of the fluid that causes 
the dust to concentrate will eventually drift past. The local dust will then 
encounter gas with a higher than normal degree of vorticity and this will 
cause the the dust concentration to begin to diminish. 

Finally, Figures [5ja-b) show the temperature and stream function in 
the full 2D domain. 



(a) Time Series: with dissipation 
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Figure 6. Representative profiles for Model C showing a prograde shear 
profile: T2 = 5.30 and <j> = 10, $ = — 10 along with b\ = —10. These 
results were generated with B = 1. (a) Time series of the absolute ampli- 
tude of F x and A. The amplitudes reach a steady average value but each 
profile experiences a small amplitude oscillation that neither grows nor de- 
cays with time. Inset is shown a higher resolution profile of the transient 
phase for the vorticity near s ~ 2. The profile goes from nominally two 
large scale humps to one. Panels (b-d) show profiles of various quantities at 
three different times: (b), initial growth phase, time s ~ 0.84, (c), transient 
adjustment phase, time s ~ 1.72, (d), steady state phase, time s ~ 20.0. In 
(d) it is clear that the vorticity takes on significant and complicated structure 
and it is also apparent that those regions where the vorticity is negative are 
strongly concentrated. 

3.3.3 Model C: Prograde Shear 

Shear does not change the basic features of the results found for Model B. 
However, it does alter the appearence of the finally developed convection 
roll by giving it a much more complicated structure. As before, the evolution 
goes through an initial growth stage followed by a transient readjustment 
phase before settling onto a roughly steady profile (Figure [SJi). The great 
difference between these results and that of Model B is that the negative 
vorticity of the roll is strongly concentrated into a small region in space 
whereas the positive vorticity, including all the substructures like in Figure 
[6jl, are spread out over larger areas of space (Figure|7p). 

ACKNOWLEDGMENTS 

I would like to thank the DPW at BRC (2002) for providing and maintaining 
the desert field facilities in which the ideas for this work first originated. 

REFERENCES 

Barge P. & Sommeria, J., 1995, A&A, 295, LI. 

Barranco J. A. & Marcus, PS., 2000, Summer Program Proceedings, Center 

For Turbulence Research, NASA Ames/Stanford Univ. 97. 
Barranco JA. & Marcus, PS., 2005, ApJ, 623, 1157. 
Bertin, G. & Cava, A., 2006, A&A, 459, 333. 

Bracco, A., Chavanis PH., Provenzale, A. & Spiegel E.A., 1999, Phys. of 

Fluids, 11, 2280. 
Chandresekhar, S., 1954, Proc. Roy. Soc. (London) A, 225, 173. 
Chapman, C. J. & Proctor, M. R. E., 1980, J. Fluid Mech., 101, 759. 
Chiang, E.I. & Goldreich, P., 1997, ApJ, 490,368. 



8 O.M. Umurhan 



(a) Temperature, s=18.9 




together these equations yield the solution 

*o = RaPafx, 



p y 1 

= T ~ 8' 



(A13) 



(b) Stream Function, s =18.9 




where Pq is designed so that there is no normal flow of fluid at the vertical 
boundaries, or in other words, *I/o = 0atj/ = ±l/2. We note that because 
of the relationships between the lowest order vorticity, stream function and 
the temperature in jA12t and IA13K the vorticity is proportional to the gra- 
dient of the temperature structure function, or in other words, Slo ~ F x . 

In order for a solution to exist for the temperature at the next order, 
a solvability condition must be enforced upon the equations at this order. 
This means that starting with, 

D 2 B 2 = -d x * -d x Oo+b 1 yd x eo+d x (<fo)De -D(* )d x e AA14) 

the solvability condition, in which the thermal flux at the vertical boundaries 
is fixed, amounts to requiring, 



Figure 7. Representative contours for Model c in the quasi-steady state 
regime, s ~ 18.9. For these contour calculations, e = 1. (a) Temperature 
fluctuation contours, 0. (b) Stream function contours, 
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D'e. 2 dy = D<3> 2 \ y y l 1 ll 2 = V- 



(A15) 



Since ©o is constant with respect to y and since *I/o is an even function of 
y, we find that the condition )A15t applied to jA14t results in a choice for 

Ro, 

rl/2 

1 + Ro / Pody = 0. (A16) 
J -1/2 

From )A13t . it follows that Ro = 12. The solution of ©2 becomes 

e 2 = T 2 f xx + M 2 (f x ) 2 + &iAT 2 / x , (A17) 

where the individual functions, T 2 , M 2 , N 2 are given by 
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D 2 T 2 = 


-1-floPo, 


T 2 = 


4 


2 ' 


(Ai8) 
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D 2 M 2 


= -RqDP , 


M 2 = 


3y 

2 


-2y\ 


(A19) 
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D 2 N 2 = z, 


N 2 = 


(A20) 



APPENDIX A: ASYMPTOTIC EXPANSIONS 



At infinite Prandtl numbers, the full governing equations A43I481 are repro- 
duced here in slightly different form and with a little bit more detail, 



(Al) 



R (l + e 2 ^jd x e~(l + e 2 8)n = e 2 f3d x V 

e 4 d T e + t 2 j(e,^) + tb 1 zd x e = e 2 d x * + t 2 d x e- 

Q = (d 2 + e 2 d 2 x ) * (A3) 

e 6 C 6 = t 6 <f>n + --- (A4) 

tt d = t 11 -^ (A5) 



these are also designed so their y-derivatives are zero at the boundaries 
y = ±1/2. Because and are order e 7 and e 6 respectively, and 
furthermore, because of the relationship between the particle divergence 
and the gas vorticity in IA41 . the lowest order expression of fA61 is 

d T 8 2 + 0Q O = 0. (A21) 

But because of the solutions jA12l and jA13l we may write 

h = 6 2 {X,T), (A22) 

in other words, the structure function <5 for the particle concentration is in- 
dependent of y at this order. The governing evolution equation for 8 2 is 
rewritten to reflect the solution Qq , 



d T 5 2 = -4>Rof x . 



(A23) 



e b d T 8 + e b {l + e 2 8)C(i = 
The following expansions are presumed 

6 = 6 + e 2 02 + e 4 04 + ■ ■ ■ (A7) 
$ = * + e 2 *2 H (A8) 

n = n + e 2 n 2 h — (A9) 

5 = 52 + ■ ■ ■ (A10) 
and similarly for the other quantities. At the lowest order )A21 is, 



D 2 B Q = 



(All) 



We may now proceed to order e 2 and solve for the next order stream 
function ty 2 . The governing equation becomes, 

D 2 * 2 = Rod x e 2 + R 2 d x e - d x <K - 8 2 U - f3d x <f (A24) 

Utilizing the solutions from the previous orders we find that the ty 2 may be 
written as 

*2 = S 2 f xxx + (R 2 - 8 2 R )P f x + Q 2 d x (f x ) 2 

, r (odd) r (even) \ £ / a nc\ 

b\L\ '-pL\ ' j fxx (A25) 
where the structure functions appearing above are given by 



which admits the solution ©o = f(X, T). Next solve the next order mo- 
mentum equation 4AH and stream function relationship {A3) 



Rod x O -n = 0, n =D*o, 



(A12) 



D 2 S 2 = R T 2 -R P , S 2 



D 2 Q 2 =R Q M 2 , Q 2 



5 4 4 160 



,(A26) 



3»/ 3 - — , (A27) 
y 40 
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D 2 L ( ° dd) =R N 2 , 
D 2 ^ 



-,2 ^(even) 



RoPq, 
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y 3 
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(A28) 



The structure functions are designed so that they are zero at the boundaries 
y = ±1/2. Next we turn to the order e 4 heat equation which is, 

d T e + J(^2,e ) + J(* , e 2 ) + b iy d x e 2 = 

8x^2 + 8 X B 2 + D 2 B 4 - (A30) 
Like at the previous order the solvability condition is, 

• 1/2 

J D 2 04<iy = J De4|^: 1 _ / 1 2 /2 = O. (A31) 

-1/2 

The solvability condition produces an evolution equation for /, namely, 

It + K fxxxx — ^(fx)x 

+ [(ra - h)fx] x + biS(fx)x + PBf X xx = 0, (A32) 
where the constants and parameters are defined by, 

.1/2 

(Sa +T 2 )dy= — , (A33) 
1/2 21 
1/2 

/>'„ / (P DM 2 )dy = - 
1/2 5 
1/2 M/2 
, j -i? 2 / Pdy + b 2 / zN 2 dy=^-^, (A35) 

1/2 J-l/2 
1/2 

S Ho / P DN 2 dy = — , (A36) 

-1/2 10 
1/2 



(A34) 



_R2 _ 

H 120 ' 



Ho 



4 even) dy: 



1/2 



1 

To' 



(A37) 



5s as the time increment between time step n and n + 1 the discretized 
equations are, 



k 



pn-1 



2<5s 
- AI 



<5.s 



MI 



(B7) 
(B8) 

(B9) 
(BIO) 

In the first discretization above, the coefficient terms in front of F^~ l and 
Af k look like Pade approximants to exponentials of C k . This is then as- 
sumed and the following replacements, 



Rearranging the terms and isolating the n + 1 time step reveals, 

25s 



? n + l 



.n+1 



l + SsC k F „_ 1 
l-5sC k k 
A™" 1 + 25sM 7 l. 



1 - SsC, 



1 + SsC k 
1 — 8sC k 



■ exp [2SsC k ] 



1 



1 - 5s£ fe 



• exp [5s£fe] 



(Bll) 



are made in (B9J- The resulting discretized equations that are simultane- 
ously solved at each time step are )B10t along with, 

= exp [25sC k ] F™- 1 + 2<5scxp [8sC k ]N£. (B12) 

Time centered differenced schemes like the simple one used here are known 
to sometimes incur a period 25s long-time numerical instability. In order to 
wash this effect out, at every 50th timestep dB12t and flBlO) are advanced 
through for one time step using an Euler type evolver as in the following, 

= wp[8sC k ]F% + 6s<sxp\5sC h \M%, (B13) 

A™ +1 = A"+<SsM£. (B14) 

This Euler advancer is also used as the very first time step of all simulations. 



APPENDIX B: NUMERICAL METHOD 

The general form of the equations we seek to solve are, 
8F 

— = CF + Af(F, A) (Bl) 
as 

8A 

— = M(F), (B2) 
ds 

where C is some linear operator and where Af is a nonlinear operator in- 
volving some combination of the arguments and AA is linear in F. 

These one dimensional equations will be solved in terms of Fourier 
expansions of F and A as in, 

F = ^F fc (s)expifcx + c.c. (B3) 

A = Afc(s) exp ik\ + c.c, (B4) 

such that each Fourier component evolves according to 

—± = C k F k +Af k (F,A) (B5) 
as 

8A k 
8s 



M k (F). (B6) 



The scalars C k are the linear operator C acting on the appropriate Fourier 
wave k and where Af k and AA k are the fcth components of each respective 
operator Af and AA . 

The symbol F n represents the nth time iterate of F. This same con- 
vention carries over to A, A4 and Af . The evolution is carried out us- 
ing a second order time centered difference scheme together with a Crank- 
Nicholson scheme for the spatial wavenumbers of the linear operator terms. 
All Afk and AA k expressions are evaluated at the centered time point. With 



